function [v] = solve_valfunc(c,s,param,glob,options)

c2           = c(glob.Nb+1:end,:);
c2           = reshape(c2, glob.n(1), glob.n(2));
v2           = griddedInterpolant(glob.B, glob.S,c2, 'linear');   

    % solve value function
B                   = menufun('bounds',s,[],param,glob,options);
obj                 = @(y)valfunc_I(v2, s,y,param,glob,options);     
y                   = goldenx(obj,B(:,1),B(:,2));

[v1]           = valfunc_I(v2,s,y,param,glob,options);

% Report Jacobian if necessary
% THIS SEEMS TO BE INCORRECT!
v2  = glob.PnK*v1;

% need to assemble their policy function p
%p     = p.*I + s(:,1).*I;

p         = menufun('p',s,y,param,glob,options);
F         = menufun('F',s,y,param,glob,options);
L        = menufun('L',s,y,param,glob,options);

adj_cost     = menufun('C',s,y,param,glob,options);

% Package output
v.v1        = v1;
v.v2        = v2;
v.p         = p;
v.y         = y;
v.F         = F;
v.l         = L;
v.adj_cost  = adj_cost;
